Данные

life_ex <- readRDS('life_expectancy_data.RDS')
str(life_ex)
## Classes 'data.table' and 'data.frame':   195 obs. of  23 variables:
##  $ Country                                : chr  "Afghanistan" "Albania" "Algeria" "Angola" ...
##  $ Year                                   : int  2019 2019 2019 2019 2019 2019 2019 2019 2019 2019 ...
##  $ Gender                                 : chr  "Female" "Female" "Female" "Female" ...
##  $ Life expectancy                        : num  66.4 80.2 78.1 64 78.1 ...
##  $ Unemployment                           : num  14.06 11.32 18.63 7.84 8.26 ...
##  $ Infant Mortality                       : num  42.9 7.7 18.6 44.5 5.1 ...
##  $ GDP                                    : num  1.88e+10 1.54e+10 1.72e+11 8.94e+10 1.69e+09 ...
##  $ GNI                                    : num  1.91e+10 1.52e+10 1.68e+11 8.19e+10 1.58e+09 ...
##  $ Clean fuels and cooking technologies   : num  36 80.7 99.3 49.6 100 ...
##  $ Per Capita                             : num  494 5396 3990 2810 17377 ...
##  $ Mortality caused by road traffic injury: num  15.9 11.7 20.9 26.1 0 ...
##  $ Tuberculosis Incidence                 : num  189 16 61 351 0 29 26 2.2 6.9 6 ...
##  $ DPT Immunization                       : num  66 99 91 57 95 ...
##  $ HepB3 Immunization                     : num  66 99 91 53 99 ...
##  $ Measles Immunization                   : num  64 95 80 51 93 ...
##  $ Hospital beds                          : num  0.432 3.052 1.8 0.8 2.581 ...
##  $ Basic sanitation services              : num  49 99.2 86.1 51.4 85.5 ...
##  $ Tuberculosis treatment                 : num  91 88 86 69 72.3 ...
##  $ Urban population                       : num  25.8 61.2 73.2 66.2 24.5 ...
##  $ Rural population                       : num  74.2 38.8 26.8 33.8 75.5 ...
##  $ Non-communicable Mortality             : num  36.2 6 12.8 19.4 17.6 ...
##  $ Sucide Rate                            : num  3.6 2.7 1.8 2.3 0.8 ...
##  $ continent                              : Factor w/ 5 levels "Africa","Americas",..: 3 4 1 1 2 2 4 2 5 4 ...
##  - attr(*, ".internal.selfref")=<externalptr> 
##  - attr(*, "sorted")= chr "Country"

Plotly график количества городского населения по странам

urb_viz <- 
  life_ex %>% 
  ggplot(aes(x = Country, y = `Urban population`, fill  = continent)) +
  facet_wrap(~continent)+
  geom_col() +
  theme(axis.text.x =element_blank())
ggplotly(urb_viz)

Plotly график количества деревенского населения по странам

rur_viz <- 
  life_ex %>% 
  ggplot(aes(x = Country, y = `Rural population`, fill  = continent)) +
  facet_wrap(~continent)+
  geom_col() +
  theme(axis.text.x = element_blank())
ggplotly(rur_viz)

T-тест на сравнение распределений Life expectancy между группами стран Африки и Америки.

ggqqplot(life_ex, x = "Life expectancy", facet.by = "continent")

stat.test <- 
  life_ex %>% 
  filter(continent == c("Africa", "Americas")) %>% 
  rstatix::t_test(`Life expectancy` ~ continent) %>% 
  add_xy_position(x = "continent")
stat.test
## # A tibble: 1 × 12
##   .y.        group1 group2    n1    n2 statistic    df       p y.position groups
##   <chr>      <chr>  <chr>  <int> <int>     <dbl> <dbl>   <dbl>      <dbl> <name>
## 1 Life expe… Africa Ameri…    26    24     -7.01  41.1 1.57e-8       86.6 <chr> 
## # ℹ 2 more variables: xmin <dbl>, xmax <dbl>
life_ex %>% 
  filter(continent == c("Africa", "Americas")) %>% 
  ggboxplot(
  y = "Life expectancy",
  x = "continent",
  ylab = "Ожидаемая продолжительность жизни", 
  xlab = "Континент", 
  add = "jitter"
  ) + 
  labs(subtitle = get_test_label(stat.test, detailed = TRUE)) + 
  stat_pvalue_manual(stat.test, tip.length = 0) 

Корреляционный анализ

new_life_ex <- 
  life_ex %>% 
  select(where(is.numeric) | where(is.integer)) %>% 
  select(-Year)
new_life_ex_cor <- cor(new_life_ex)
corrplot(corr = new_life_ex_cor,
         method = "color",
         order = "hclust",
         tl.cex = 0.5)

ggcorr(new_life_ex, 
       label = T, 
       label_size = 2, 
       label_round = 1, 
       geom = "circle",
       size = 2,
       hjust = 0.5,
       angle = -30)

Heatmap и иерархическая кластеризация

Мной было предположено, что количество кластеров может быть 5 по числу рассматриваемых континентов, поскольку это во многом определяет уровень медицины, экокномическое развитие, климатические условия и так далее. Кластеризация была проведена по строкам и по колонкам, для каждого из 5 кластеров можно сделать вывод какие переменные отличают его среди других кластеров, а также можно сделать выводы о коллинеарности части переменных, подобно корреляции. Первый кластер включает в себя небольшое количество наблюдений и для него характерны наиболее низкие значения по Measles, DPT и HepB3 Immunization. Для второго кластера, по сравнению со всеми остальными, характерны высокие значения Mortality caused by road traffic injury, Tuberculosis Incidence, Infant Mortality. Третий кластер самый малочисленный, но выделяется высокими значениями скоррелированных GDP и GNI. Четвертый кластер имеет наиболее низкие значения Mortality caused by road traffic injury, Tuberculosis Incidence, Infant Mortality, но высокие Clean fuels and cooking technologies, Per Capita, Urban population.

new_life_ex_scaled <- scale(new_life_ex)
pheatmap(new_life_ex_scaled, 
         show_rownames = FALSE, 
         clustering_distance_rows = "euclidean",
         clustering_method = "ward.D2", 
         cutree_rows = 5,
         cutree_cols = length(colnames(new_life_ex_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

PCA анализ

Проводим метод главных компонент и визуализируем Cumulative Proportion, видим, что первая главная компонента отвечает за наибольшую дисперсию в наблюдениях. Выбираем первые две компаненты для визуализации и смотрим с какими переменными свазаны наши главные оси. С точки зрения вариации наиболее значимы переменные Life expectancy и Measles, DPT и HepB3 Immunization, причем Life expectancy сильно связана с PC1, что видно по углу наклона, то есть наблюдения с высоким Life expectancy будут лежать в положительной части PC1. Measles, DPT и HepB3 Immunization связаны как с PC1, так и с PC2, при этом Measles, DPT и HepB3 Immunization прямо пропорциональны PC1 и обратно пропорциональны PC2, то есть наблюдения с высокми значениями по Measles, DPT и HepB3 Immunization будут лежать в 4 четверти координатной плоскости. Если смотреть на то, какие переменные вносят наибольший вклад в PC1 и PC2, то в случае PC1 это Life expectancy, Infant Mortality, Basic sanitation services, Clean fuels and cooking technologies, в случае PC2 это Measles, DPT и HepB3 Immunization.

new_life_ex.pca <- prcomp(new_life_ex, 
                          scale = T)

summary(new_life_ex.pca)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6    PC7
## Standard deviation     2.7526 1.4841 1.3952 1.17177 1.08375 0.96347 0.9288
## Proportion of Variance 0.3988 0.1159 0.1025 0.07227 0.06182 0.04886 0.0454
## Cumulative Proportion  0.3988 0.5147 0.6172 0.68945 0.75126 0.80012 0.8455
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.85740 0.69263 0.68937 0.59106 0.54986 0.47085 0.36596
## Proportion of Variance 0.03869 0.02525 0.02501 0.01839 0.01591 0.01167 0.00705
## Cumulative Proportion  0.88421 0.90946 0.93447 0.95286 0.96877 0.98044 0.98749
##                           PC15    PC16    PC17    PC18      PC19
## Standard deviation     0.34546 0.26941 0.20224 0.06968 7.273e-16
## Proportion of Variance 0.00628 0.00382 0.00215 0.00026 0.000e+00
## Cumulative Proportion  0.99377 0.99759 0.99974 1.00000 1.000e+00
fviz_eig(new_life_ex.pca, addlabels = T)

fviz_pca_var(new_life_ex.pca, col.var = "contrib")

fviz_pca_var(new_life_ex.pca, 
             select.var = list(contrib = 4),
             col.var = "contrib")

fviz_contrib(new_life_ex.pca, choice = "var", axes = 1, top = 24)

fviz_contrib(new_life_ex.pca, choice = "var", axes = 2, top = 24)

PCA визуализация

На данном биплоте мы можем проанализировать где наши наблюдения находятся в новой координатной плоскости, а также понять какие для них характерны значения по тем или иным переменным, основываясь на том, как переменные соотносятся с главными компонентами. Проведя визуализацию наблюдений цветом по континентам, мы в целом можем сделать выводы насколько различаются и различаются ли вообще наблюдения по странам с одного континента. В данном случаем мы видим плотно сгруппированные кластеры наблюдений по Европе и Америке, которые при этом не сильно различаются между собой, а также кластер Африканских стран, который размазан больше, но при этом достаточно отделяется от других кластеров. Наблюдения из Азии и Океании плохо разделяются на кластеры и отдельные наблюдения достаточно сильно отдалены друг от друга, что говорит о то, что часть стран по разным переменным схожа с Африканскими странами, а часть с Европой/Америкой.

pca_life_ex <- as.data.frame(new_life_ex.pca$x)

bip_life_ex <- 
  ggbiplot(new_life_ex.pca,
           scale=F,
           labels = F,
           groups = as.factor(life_ex$continent), 
           ellipse = T,
           alpha = 0,
           varname.abbrev = T)
bip_life_ex <- 
  bip_life_ex + 
  geom_point(data = pca_life_ex, 
                    aes(x = PC1, 
                        y = PC2, 
                        color = life_ex$continent, 
                        text = life_ex$Country, alpha = 0.5)) +
  labs(color = "Continent")

ggplotly(bip_life_ex, tooltip = c("text"))

UMAP

Несмотря на другое располжение наблюдений на кординатной плоскости методом UMAP относительно PCA, тенденция разделения наблюдений на группы по континентам схожая, хорошо выделяется кластер наблюдений стран Европа/Америка, а также дасточно сильно удаленный от них кластер стран Африки, страны Океании и Азии попадают как в Африканский кластер, так и кластер Европа/Америка, что объясняется тем, что условия в данных странах могут быть приближены, как к странам Африки, так и к странам Европы/Америки, что во многом связано и с многочисленностью стран в данных регионах и значительной географической протяженностью данных континентов.

umap_prep <- recipe(~., data = new_life_ex) %>% 
  step_normalize(all_predictors()) %>%
  step_umap(all_predictors()) %>%
  prep() %>%
  juice()

umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) +
  geom_point(aes(color = as.factor(life_ex$continent)),
             alpha = 0.7, size = 2) +
  labs(color = NULL)